home *** CD-ROM | disk | FTP | other *** search
/ Monster Media 1996 #15 / Monster Media Number 15 (Monster Media)(July 1996).ISO / prog_c / cuj0696.zip / DWYER.ZIP / LIB / POWI.C < prev    next >
C/C++ Source or Header  |  1995-12-12  |  3KB  |  163 lines
  1. /*                            powi.c
  2.  *
  3.  *    Real raised to integer power
  4.  *
  5.  *
  6.  *
  7.  * SYNOPSIS:
  8.  *
  9.  * double x, y, powi();
  10.  * int n;
  11.  *
  12.  * y = powi( x, n );
  13.  *
  14.  *
  15.  *
  16.  * DESCRIPTION:
  17.  *
  18.  * Returns argument x raised to the nth power.
  19.  * The routine efficiently decomposes n as a sum of powers of
  20.  * two. The desired power is a product of two-to-the-kth
  21.  * powers of x.  Thus to compute the 32767 power of x requires
  22.  * 28 multiplications instead of 32767 multiplications.
  23.  *
  24.  *
  25.  *
  26.  * ACCURACY:
  27.  *
  28.  *
  29.  *                      Relative error:
  30.  * arithmetic   x domain   n domain  # trials      peak         rms
  31.  *    DEC       .04,26     -26,26    100000       2.7e-16     4.3e-17
  32.  *    IEEE      .04,26     -26,26     50000       2.0e-15     3.8e-16
  33.  *    IEEE        1,2    -1022,1023   50000       8.6e-14     1.6e-14
  34.  *
  35.  * Returns MAXNUM on overflow, zero on underflow.
  36.  *
  37.  */
  38.  
  39. /*                            powi.c    */
  40.  
  41. /*
  42. Cephes Math Library Release 2.1:  January, 1989
  43. Copyright 1984, 1987, 1989 by Stephen L. Moshier
  44. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  45. */
  46.  
  47. #include "mconf.h"
  48. extern double MAXNUM, MAXLOG, MINLOG, LOGE2;
  49.  
  50. double powi( x, nn )
  51. double x;
  52. int nn;
  53. {
  54. int n, e, sign, asign, lx;
  55. double w, y, s;
  56. double log(), frexp();
  57.  
  58. if( x == 0.0 )
  59.     {
  60.     if( nn == 0 )
  61.         return( 1.0 );
  62.     else if( nn < 0 )
  63.         return( MAXNUM );
  64.     else
  65.         return( 0.0 );
  66.     }
  67.  
  68. if( nn == 0 )
  69.     return( 1.0 );
  70.  
  71.  
  72. if( x < 0.0 )
  73.     {
  74.     asign = -1;
  75.     x = -x;
  76.     }
  77. else
  78.     asign = 0;
  79.  
  80.  
  81. if( nn < 0 )
  82.     {
  83.     sign = -1;
  84.     n = -nn;
  85.     }
  86. else
  87.     {
  88.     sign = 1;
  89.     n = nn;
  90.     }
  91.  
  92. /* Overflow detection */
  93.  
  94. /* Calculate approximate logarithm of answer */
  95. s = frexp( x, &lx );
  96. e = (lx - 1)*n;
  97. if( (e == 0) || (e > 64) || (e < -64) )
  98.     {
  99.     s = (s - 7.0710678118654752e-1) / (s +  7.0710678118654752e-1);
  100.     s = (2.9142135623730950 * s - 0.5 + lx) * nn * LOGE2;
  101.     }
  102. else
  103.     {
  104.     s = LOGE2 * e;
  105.     }
  106.  
  107. if( s > MAXLOG )
  108.     {
  109.     mtherr( "powi", OVERFLOW );
  110.     y = MAXNUM;
  111.     goto done;
  112.     }
  113.  
  114. #if 0
  115. if( s < MINLOG )
  116.     return(0.0);
  117.  
  118. /* Handle tiny denormal answer, but with less accuracy
  119.  * since roundoff error in 1.0/x will be amplified.
  120.  * The precise demarcation should be the gradual underflow threshold.
  121.  */
  122. if( (s < (-MAXLOG+2.0)) && (sign < 0) )
  123.     {
  124.     x = 1.0/x;
  125.     sign = -sign;
  126.     }
  127. #else
  128. /* do not produce denormal answer */
  129. if( s < -MAXLOG )
  130.     return(0.0);
  131. #endif
  132.  
  133.  
  134. /* First bit of the power */
  135. if( n & 1 )
  136.     y = x;
  137.  
  138. else
  139.     {
  140.     y = 1.0;
  141.     asign = 0;
  142.     }
  143.  
  144. w = x;
  145. n >>= 1;
  146. while( n )
  147.     {
  148.     w = w * w;    /* arg to the 2-to-the-kth power */
  149.     if( n & 1 )    /* if that bit is set, then include in product */
  150.         y *= w;
  151.     n >>= 1;
  152.     }
  153.  
  154.  
  155. done:
  156.  
  157. if( asign )
  158.     y = -y; /* odd power of negative number */
  159. if( sign < 0 )
  160.     y = 1.0/y;
  161. return(y);
  162. }
  163.